En este notebook, te mostraremos cómo extraer características en modo batch utilizando el enfoque de la interfaz de línea de comandos (cli).
Utilizaremos datos del repositorio https://github.com/rcuocolo/PROSTATEx_masks
Primero, necesitamos verificar e instalar las bibliotecas necesarias para realizar estos pasos más fácilmente.
!pip install --progress-bar off -q gdown!pip install --upgrade --pre SimpleITK --find-links https://github.com/SimpleITK/SimpleITK/releases/tag/latest --progress-bar off -q!pip install --progress-bar off -q matplotlib!pip install --progress-bar off -q pyradiomics!pip install --progress-bar off -q pyyaml!pip install --progress-bar off -q pandas
Installing build dependencies ... done
Getting requirements to build wheel ... done
Preparing metadata (pyproject.toml) ... done
Preparing metadata (setup.py) ... done
Preparing metadata (setup.py) ... done
Building wheel for pyradiomics (setup.py) ... done
Building wheel for docopt (setup.py) ... done
Importar librerías necesarias
import osimport pandas as pdimport csvimport yamlimport SimpleITK as sitkimport radiomicsimport matplotlib.pyplot as pltimport gdownimport randomimport numpy as npimport ipywidgets as widgetsfrom ipywidgets import interact, fixed, Layout, Button, HBoxfrom IPython.display import display, clear_outputimport seaborn as snsfrom sklearn.feature_selection import VarianceThreshold
# prompt: Montar mi unidad de google drivefrom google.colab import drivedrive.mount('/content/drive')
Mounted at /content/drive
Descargar imágenes ProstateX y archivo de extracción de máscaras y parámetros
import osifnot os.path.isdir('input_t2w.csv'):!gdown 1dNqdVsN0PlCG01sl6wUCh9QJ5NnkOrDS-q!unzip lesions_T2.zip>/dev/null!rm lesions_T2.zip>/dev/null!rm -r __MACOSX >/dev/nullelse:print("lesions_T2 folder already exists, the zip was not downloaded")ifnot os.path.isfile('Params_T2_BW.yaml'):!gdown 12VOQCPFK4HARrwQzbB4fr_f2PepGPFWm-qelse:print("Params_T2_BW.yaml already exists, the file was not downloaded")print('Done!')
Done!
Visualizar una imagen al azar y la máscara de la lesión
# Rutas a las carpetas (relativas para Colab)ruta_imagenes ='/content/lesions_T2/Images'# Usamos rutas relativas para Colabruta_mascaras ='/content/lesions_T2/Masks'# Usamos rutas relativas para colab# Variables globalespareja_actual =Noneparejas_archivos = []output_widget = widgets.Output() # Widget de salida para mostrar las imágenesdef mostrar_imagen_y_mascara(ruta_imagen, ruta_mascara, zoom_x=0, zoom_y=0, zoom_size=50):""" Muestra las tres imágenes (MRI, MRI+máscara, MRI+máscara con zoom). """try: imagen_sitk = sitk.ReadImage(ruta_imagen) mascara_sitk = sitk.ReadImage(ruta_mascara)exceptRuntimeErroras e:print(f"Error al leer: {e}")return imagen_array = sitk.GetArrayFromImage(imagen_sitk) mascara_array = sitk.GetArrayFromImage(mascara_sitk)if imagen_array.shape != mascara_array.shape:iflen(imagen_array.shape) ==3andlen(mascara_array.shape) ==4:if mascara_array.shape[3] ==1: mascara_array = np.squeeze(mascara_array, axis=3)else:print("Dimensiones incompatibles (4D).")returneliflen(imagen_array.shape) !=3orlen(mascara_array.shape) !=3:print("Dimensiones incompatibles.")return cortes_con_lesion = np.where(np.any(mascara_array !=0, axis=(1, 2)))[0]if cortes_con_lesion.size ==0:print("No se encontraron cortes con lesión.")iflen(imagen_array.shape) ==3: corte = random.randint(0, imagen_array.shape[0] -1)else: imagen_slice = imagen_array mascara_slice = mascara_arrayelse: corte = random.choice(cortes_con_lesion)iflen(imagen_array.shape) ==3: imagen_slice = imagen_array[corte, :, :] mascara_slice = mascara_array[corte, :, :]# --- Preparación para el zoom --- zoom_x =max(zoom_size, min(zoom_x, imagen_slice.shape[1] - zoom_size -1)) zoom_y =max(zoom_size, min(zoom_y, imagen_slice.shape[0] - zoom_size -1)) fila_inicio =max(0, zoom_y - zoom_size) fila_fin =min(imagen_slice.shape[0], zoom_y + zoom_size) col_inicio =max(0, zoom_x - zoom_size) col_fin =min(imagen_slice.shape[1], zoom_x + zoom_size) imagen_zoom = imagen_slice[fila_inicio:fila_fin, col_inicio:col_fin] mascara_zoom = mascara_slice[fila_inicio:fila_fin, col_inicio:col_fin] mascara_zoom_rgba = np.zeros((mascara_zoom.shape[0], mascara_zoom.shape[1], 4)) mascara_zoom_rgba[mascara_zoom !=0, :3] = plt.cm.jet(mascara_zoom[mascara_zoom !=0])[:, :3] mascara_zoom_rgba[mascara_zoom !=0, 3] =0.5# --- Creación de la figura y subplots DENTRO del contexto del widget de salida ---with output_widget: clear_output(wait=True) # Borrar la salida anterior fig, axes = plt.subplots(1, 3, figsize=(18, 6))# --- Mostrar imagen original (sin zoom) --- axes[0].imshow(imagen_slice, cmap='gray') axes[0].set_title(f'Imagen MRI (Corte: {corte})') axes[0].axis('off')# --- Mostrar imagen con máscara superpuesta (sin zoom) --- mascara_rgba = np.zeros((mascara_slice.shape[0], mascara_slice.shape[1], 4)) mascara_rgba[mascara_slice !=0, :3] = plt.cm.jet(mascara_slice[mascara_slice !=0])[:, :3] mascara_rgba[mascara_slice !=0, 3] =0.5 axes[1].imshow(imagen_slice, cmap='gray') axes[1].imshow(mascara_rgba) axes[1].set_title(f'Máscara Superpuesta (Corte: {corte})') axes[1].axis('off')# --- Mostrar imagen con máscara superpuesta (zoom) --- axes[2].imshow(imagen_zoom, cmap='gray') axes[2].imshow(mascara_zoom_rgba) axes[2].set_title('Máscara Superpuesta (Zoom)') axes[2].axis('off') plt.tight_layout() plt.show()def obtener_parejas_archivos(ruta_imagenes, ruta_mascaras):"""Obtiene y empareja archivos.""" imagenes =sorted([f for f in os.listdir(ruta_imagenes) if f.endswith('.nii.gz') and"ProstateX"in f]) mascaras =sorted([f for f in os.listdir(ruta_mascaras) if f.endswith('.nii.gz') and"ProstateX"in f]) parejas = []for img in imagenes: img_id = img.split('_')[0]for masc in mascaras:if img_id in masc: ruta_img = os.path.join(ruta_imagenes, img) ruta_masc = os.path.join(ruta_mascaras, masc) parejas.append((ruta_img, ruta_masc))break# Importante: salir del bucle interno una vez encontrada la parejareturn parejasdef actualizar_imagen(b):"""Función para actualizar la imagen al hacer clic en el botón."""global pareja_actual, parejas_archivosifnot parejas_archivos: parejas_archivos = obtener_parejas_archivos(ruta_imagenes, ruta_mascaras) pareja_actual = random.choice(parejas_archivos) imagen_sitk = sitk.ReadImage(pareja_actual[0]) imagen_array = sitk.GetArrayFromImage(imagen_sitk)iflen(imagen_array.shape) ==3: dim_x = imagen_array.shape[2] dim_y = imagen_array.shape[1]else: dim_x = imagen_array.shape[1] dim_y = imagen_array.shape[0] zoom_x_slider.max= dim_x -1 zoom_x_slider.value = dim_x //2 zoom_y_slider.max= dim_y -1 zoom_y_slider.value = dim_y //2 zoom_size_slider.max=min(dim_x, dim_y) //2 zoom_size_slider.value =50 mostrar_imagen_y_mascara(pareja_actual[0], pareja_actual[1], zoom_x_slider.value, zoom_y_slider.value, zoom_size_slider.value)# --- Configuración inicial ---parejas_archivos = obtener_parejas_archivos(ruta_imagenes, ruta_mascaras)ifnot parejas_archivos:print("No se encontraron parejas de archivos.")#Si no se encuentran parejas, no se puede continuar. Se podría poner un mensaje en un widget.else: pareja_actual = random.choice(parejas_archivos) imagen_sitk_inicial = sitk.ReadImage(pareja_actual[0]) imagen_array_inicial = sitk.GetArrayFromImage(imagen_sitk_inicial)iflen(imagen_array_inicial.shape) ==3: dim_x = imagen_array_inicial.shape[2] dim_y = imagen_array_inicial.shape[1]else: #Si es 2D dim_x = imagen_array_inicial.shape[1] dim_y = imagen_array_inicial.shape[0]# --- Crear widgets --- zoom_x_slider = widgets.IntSlider(min=0, max=dim_x -1, step=1, value=dim_x //2, description='Zoom X', continuous_update=False, layout=Layout(width='50%')) zoom_y_slider = widgets.IntSlider(min=0, max=dim_y -1, step=1, value=dim_y //2, description='Zoom Y', continuous_update=False, layout=Layout(width='50%')) zoom_size_slider = widgets.IntSlider(min=10, max=min(dim_x,dim_y) //2, step=1, value=50, description='Zoom Size', continuous_update=False, layout=Layout(width='50%')) boton_cambiar = Button(description="Cambiar Imagen") boton_cambiar.on_click(actualizar_imagen)# --- Configurar la interfaz --- ui = HBox([boton_cambiar, zoom_x_slider, zoom_y_slider, zoom_size_slider])# Ya no usamos interactive_output, sino el widget de salida directamente display(ui, output_widget)# Mostrar la imagen inicial mostrar_imagen_y_mascara(pareja_actual[0], pareja_actual[1], zoom_x_slider.value, zoom_y_slider.value, zoom_size_slider.value)
Veamos nuestro archivo de parámetros
Seleccione el icono de la carpeta en la barra de la izquierda y haga doble clic en Params_T2_BW.yaml para verlo.
Normalización
Dado que estamos trabajando con imágenes T2w MRI vamos a normalizar las imágenes utilizando el parámetro normalize que establecemos como true. Esto fija la intensidad media de las imágenes en cero y escala la desviación estándar para que sea igual a uno. A continuación escalamos la desviación estándar de las intensidades de las imágenes a 100 estableciendo normalizeScale en true.
El cálculo de algunas características (a saber, la energía, la energía total y RMS) se ve afectada por la existencia de valor negativo y ya que hemos normalizado nuestras imágenes para tener media 0 y desviación estándar de 100 tenemos que asegurarnos de que la mayoría de nuestras intensidades de imagen normalizada es igual o superior a 0. Podemos hacer esto usando el voxelArrayShift. Seleccionamos un valor de 300, que asumiendo una distribución normal de las intensidades de la imagen aseguraría que el 99.7% de ellas están por encima de 0.
Ahorro de memoria
El valor preCrop ajustado a true nos permite ahorrar algo de memoria durante el cálculo de las características (y posiblemente evitar errores).
Tipo de extracción de características (2D vs 3D)
Tenemos un conjunto de datos en el que nuestras imágenes son anisótropas (las dimensiones de los vóxeles no son las mismas para los tres ejes) y, más concretamente, el grosor del corte es mucho mayor que la resolución en el plano. Debido a esto, realizaremos la extracción de características radiómicas 2D en lugar de 3D. Ya que nuestras imágenes son Axiales T2w, necesitamos definir que deseamos realizar una extracción de características 2D estableciendo force2D a true y estableciendo el eje en el plano a 0 usando force2Ddimension (Axial - 0, Coronal - 1, Sagital - 2).
Interpolación/Remuestreo
Para asegurar que estamos extrayendo características radiómicas que son comparables entre diferentes imágenes necesitamos estandarizar las dimensiones de los vóxeles en el plano. Hacemos esto ajustando el interpolator a sitkBSpline y el resampledPixelSpacing a una resolución elegida [0.6, 0.6, 0] (el 0 en la tercera posición significa que no remuestrearemos el grosor del corte y sólo la resolución en el plano de las imágenes será transformada a 0.6 mm x 0.6 mm, lo cual está bien porque las características serán extraídas en 2D).
Resegmentación
La resegmentación permite redefinir la máscara para, por ejemplo, descartar intensidades en nuestra máscara que estén fuera de \(mean\ \pm\ 3\times std\). Podemos hacerlo ajustando resegmentMode a sigma y resegmentRange a [-3, 3].
Discritización de la intensidad
Para asegurarnos de que el número de bins utilizados para la extracción de características está entre 30 y 130, realizamos una primera extracción de características para obtener el original_firstorder_Range de nuestros pacientes y seleccionamos el binWidth para que sea mayor que 30 y menor que 130 en todos los pacientes.
Filtros y características
Utilizamos imageType para definir los tipos de imágenes de las que queremos extraer las características radiómicas y featureClass para seleccionar el tipo de características a extraer.
Realizar la extracción de características radiómicas - Crear el fichero de entrada
Para la extracción de características radiómicas utilizaremos el modo batch de PyRadiomics, donde tendremos que proporcionar como entrada el fichero .csv que contiene 3 columnas. La primera columna es el ID, la segunda la Image y la tercera la Mask, como se muestra a continuación.
ifnot os.path.isfile('input_t2w.csv'): dir_image_path ='lesions_T2/Images' dir_mask_path ='lesions_T2/Masks' img_files =sorted([os.path.join(dir_image_path, x) for x in os.listdir(dir_image_path) if'.nii'in x]) msk_files =sorted([os.path.join(dir_mask_path, x) for x in os.listdir(dir_mask_path) if'.nii.gz'in x]) input_csv_t2 = [["ID", "Image", "Mask"]]for msk_file in msk_files: patient_id = os.path.basename(msk_file)[:14] image = [x for x in img_files if patient_id in x][0] input_csv_t2.append([patient_id, image, msk_file])withopen('input_t2w.csv', 'w', newline='') as f: writer = csv.writer(f) writer.writerows(input_csv_t2)
# prompt: crear una tabla simplificada de input_df tomando solamente 20 registros al azar y guardarla como csv con el mismo formato que input_t2w.csv# Lee el archivo CSV de entradainput_df = pd.read_csv('input_t2w.csv')# Muestra las primeras 5 filas del DataFrame para verificarprint(input_df.head())# Selecciona 2 filas aleatoriassimplified_df = input_df.sample(n=2)# Guarda el DataFrame simplificado como un nuevo archivo CSVsimplified_df.to_csv('simplified_input.csv', index=False) # index=False evita guardar el índiceprint("Archivo 'simplified_input.csv' creado con éxito.")
Entradas necesarias para la ejecución por lotes de PyRadiomics
Además del fichero de entrada tendremos que proporcionar la ruta del fichero de salida donde se guardarán los valores de las características de las diferentes imágenes (usar -o antes de indicar la ruta), definir que queremos la salida como un fichero csv añadiendo -f csv, y finalmente indicar la ruta del fichero de parámetros usando -p.
Ejecutar la extracción por lotes de características radiómicas
Se aconseja ejecutar la siguiente línea en su terminal o línea de comandos para acelerar la extracción. Para ello, elimine el «!» inicial de la siguiente línea
[2025-02-14 20:57:16] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:16] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:16] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:17] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:17] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:18] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:18] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:19] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:19] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:19] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:21] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:21] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:23] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:24] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:24] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:28] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:29] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:30] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:31] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:32] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:35] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:36] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:40] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:41] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:57:42] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:58:10] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:58:10] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:58:21] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:58:35] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
[2025-02-14 20:58:35] W: radiomics.glcm: GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated
Comprobemos el fichero de salida radiomics
Como se puede ver, contiene las mismas tres columnas (ID, Imagen, Máscara) de nuestra entrada, diez columnas que proporcionan información sobre las versiones, las propiedades de la imagen y los ajustes de extracción de características, seguidas de los valores de las características radiómicas extraídas.
# prompt: Grabar el fichero en drive con la ruta output_file_pathimport shutil# ... (rest of your code)# Assuming you have already mounted your Google Driverepo_path ='/content/drive/MyDrive/radiomics/'# Replace with the correct pathoutput_file_path = os.path.join(repo_path, 'output_t2w.csv')# Use shutil to copy the file. This is more robust than relying on shell commands.try: shutil.copyfile('output_t2w.csv', output_file_path)print(f"File 'output_t2w.csv' copied successfully to {output_file_path}")exceptFileNotFoundError:print("Error: 'output_t2w.csv' not found. Ensure the radiomics extraction completed successfully.")exceptExceptionas e:print(f"An error occurred: {e}")
File 'output_t2w.csv' copied successfully to /content/drive/MyDrive/radiomics/output_t2w.csv
# prompt: De la tabla anterior, hacer un gráfico de correlación de los valores que son numéricosimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsfrom sklearn.feature_selection import VarianceThresholdradiomics_df = pd.read_csv('output_t2w.csv')for col in radiomics_df.select_dtypes(include=['object']).columns:try: radiomics_df[col] = pd.to_numeric(radiomics_df[col], errors='coerce')except (ValueError, TypeError):passnumeric_cols = radiomics_df.select_dtypes(include=['number']).dropna(axis=1, how='all')if numeric_cols.empty:print("No numeric columns")else:# --- Filtrado por varianza --- selector = VarianceThreshold(threshold=0.01) # Elimina características con varianza < 0.01. ¡Ajusta este umbral! selected_features = selector.fit_transform(numeric_cols) selected_cols = numeric_cols.columns[selector.get_support()] # Obtiene los nombres de las columnas seleccionadas filtered_df = pd.DataFrame(selected_features, columns=selected_cols) output_file_path = os.path.join(repo_path, 'coor_file.csv') filtered_df.to_csv(output_file_path, index=False)# --- Ahora calcula la correlación y grafica SOLO con las características filtradas --- correlation_matrix = filtered_df.corr() output_file_path = os.path.join(repo_path, 'correlation_matrix.csv') correlation_matrix.to_csv(output_file_path, index=False)
# prompt: cargar un archivo csv con pandasimport pandas as pd# Assuming 'input_t2w.csv' is in your current Colab environmenttry: correlation_matrix = pd.read_csv('correlation.csv')print(correlation_matrix.head()) # Display the first few rowsexceptFileNotFoundError:print("Error: 'input_t2w.csv' not found. Please upload the file or check the file path.")except pd.errors.EmptyDataError:print("Error: 'input_t2w.csv' is empty.")except pd.errors.ParserError:print("Error: Could not parse 'input_t2w.csv'. Check file format.")exceptExceptionas e:print(f"An unexpected error occurred: {e}")
# prompt: Crear un gráfico de correlación con plotlyimport plotly.express as px# Assuming 'correlation_matrix' is your DataFramefig = px.imshow(correlation_matrix, labels=dict(x="Feature", y="Feature", color="Correlation"), x=correlation_matrix.columns, y=correlation_matrix.columns, color_continuous_scale='RdBu', # Use a diverging color scale zmin=-1, zmax=1) # Set color scale limits for better visualizationfig.update_layout(title='Correlation Matrix Heatmap', width=800, height=800)fig.show()
